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Abstract 

We present a numerical scheme that can be combined with any fixed boundary finite 
element based Poisson or Grad-Shafranov solver to compute the first and second partial 
derivatives of the solution to these equations with the same order of convergence as the 
solution itself. At the heart of our scheme is an efficient and accurate computation of 
the Dirichlet to Neumann map through the evaluation of a singular volume integral and 
the solution to a Fredholm integral equation of the second kind. Our numerical method 
is particularly useful for magnetic confinement fusion simulations, since it allows the 
evaluation of quantities such as the magnetic field, the parallel current density and the 
magnetic curvature with much higher accuracy than has been previously feasible on the 
affordable coarse grids that are usually implemented. 

Keywords: plasma, equilibrium, magnetic confinement fusion, Grad-Shafranov 
equation, finite elements, integral equations, quadrature by expansion 


1. Introduction 

In computational physics, one often computes fields by expressing them in terms of 
a potential or stream function and then solving the resulting elliptic partial differential 
equation for the potential or stream function. The most common examples are problems 
involving electrostatic fields. The electric field is expressed in terms of the electric po¬ 
tential, which is then found via Poisson’s equation. Another common situation is one 
in which the Euler equations for an incompressible fluid are expressed in terms of vor- 
ticity and a stream function: the desired velocity field is obtained from derivatives of 
the stream function, which is computed by solving Poisson’s equation. A third impor¬ 
tant example, and the focus of this article, occurs in magnetic confinement fusion. The 
confining magnetic field in toroidally axisymmetric geometries is expressed in terms of 
a stream function associated with the magnetic flux, and the equilibrium magnetic con¬ 
figuration is computed by solving an elliptic PDE for this stream function known as the 
Grad-Shafranov equation nq. 

When potentials and stream functions are computed using conventional finite differ¬ 
ence or finite element schemes, numerical derivatives must be calculated to evaluate the 
fields. As a result, the convergence of the field solution is at least one order lower than 
the convergence order of the potential or stream function [anU- If derivatives of the 
fields themselves are required, at least two orders of convergence are lost. In simulations 
of magnetic fusion plasmas, this loss of accuracy is a particularly acute issue. Gertain 
physical parameters that have a key role on the stability and transport properties of hot 
plasmas, such as the parallel current density and the magnetic curvature, depend on first 
derivatives of the magnetic field or equivalently second derivatives of the flux. Popular 
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magnetic equilibrium solvers based on a finite element formulation of the Grad-Shafranov 
equation compute the magnetic stream function with a numerical error that decreases 
as N~'^, where N is the number of grid points in either direction of the two dimensional 
problem [SlElCli. The numerical errors for the derivatives of the magnetic field thus 
only decrease as N~‘^. 

There are several ways to mitigate this loss of accuracy. One can for example use 
higher order finite elements, as has been effectively implemented in |1]. However, in doing 
so, one increases the number of degrees of freedom, which leads to a larger computational 
cost. One can also use a spectral representation for the solution [9] or a combined modal- 
Green’s function representation [10]. It is indeed well known that with these methods, 
the loss of accuracy in the computation of derivatives is not an order of convergence, 
but instead a constant, as demonstrated in mm for the Grad-Shafranov equation. 
While satisfactory from a computational point of view, a weakness of such an approach 
is that the solvers rely on grids that are not typically the grids necessitated by existing 
stability, transport, or wave heating codes that take the computed equilibrium as an input 
[HI [131 [IS ESI- This is the motivation for the present work. We describe a method that 
combines the accuracy properties of a Green’s function formulation with the advantages 
of relying on a popular finite element solver for the computation of the solution to the 
Grad-Shafranov equation. The end result is the evaluation of first and second derivatives 
of the flux that has the same numerical accuracy and convergence properties as the flux 
itself while using a physically relevant grid, without having to increase the number of 
grid points. 

Our basic idea is to obtain linear partial differential equations for the first and second 
partial derivatives of the potential by analytically differentiating the original PDE (say 
the Grad-Shafranov equation or Poisson’s equation). For instance, in the case of Poisson’s 
equation in 2-D we see that the equations for the first derivatives are given by 

Au{x,y) = f{x,y) Au^ = f^{x,y), Auy = fy{x,y). (1) 

Instead of computing the derivatives of u from the output of the finite element solver, we 
use that same solver to solve the linear partial differential equations for the derivatives 
of u. 

One might then think that using the same elements as the ones used to compute 
u automatically leads to the same order of convergence for Vu as for u. This is only 
partially true in general because one also needs the boundary condition for the normal 
derivative of u on the boundary of the computational domain to the same accuracy as u 
itself. The simplest way to accomplish this task is by taking normal numerical derivatives 
of u on the boundary. This, however, leads to the loss of one order in accuracy for each 
numerical derivative taken so that we are no better off than we were by simply taking 
numerical derivatives of the original solution over the whole domain. In other words, what 
is needed on the boundary is a method to compute the Dirichlet to Neumann map that 
does not lead to any loss in the order of convergence of the solution. The development 
of such a procedure is the main new contribution of the present work. 

Specifically, we construct a numerical scheme that achieves the desired goal as follows. 
First, we re-express the Grad-Shafranov equation as a semi-linear Poisson equation with 
source function F(x, u). We then decompose the solution to Poisson’s equation as the sum 
of a particular solution that does not satisfy the proper Dirichlet boundary condition 
in general, plus a homogeneous solution that solves Laplace’s equation and is chosen 
so that the full solution satisfies the proper boundary condition. We write as the 
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volume integral of the product of -F(x, u) with the free-space Green’s function for Poisson’s 
equation in two dimensions. By differentiating under the integral sign, we have an exact 
expression for the normal derivative of uP which may be evaluated using a high-order 
quadrature. To compute the normal derivatives without ever evaluating derivatives 
normal to the boundary numerically, we introduce the harmonic conjugate U of and 
use Green’s second identity to derive a Fredholm integral equation of the second kind for 
U on the boundary of the computational domain. After solving this integral equation 
for U to the same order of accuracy as u, we can compute spectrally accurate tangential 
derivatives of U on the boundary. Since U is the conjugate gradient of this is equivalent 
to computing normal derivatives of with spectral accuracy, which is precisely what is 
needed. Note that while this work is mainly motivated by magnetic fusion applications 
and aimed at improving the accuracy of hnite element based Grad-Shafranov solvers, our 
numerical method relies on a reformulation of the Grad-Shafranov equation as a Poisson 
problem, and can therefore also be implemented in combination with Poisson solvers. 

Importantly, our approach is independent of the particular hnite element solver and 
grid used in solving the original elliptic PDF. However, the use of the harmonic conjugate 
means the approach is limited to two-dimensional problems. The method also requires 
that accurate representations of the derivatives of F’(x, u) be available, and that the 
boundary of the computational domain be smooth. These constraints are satished for a 
large number of physically relevant problems. The extension to 3-D remains a topic for 
future research. 

The structure of the article is as follows. In Section we describe in more detail 
the general philosophy guiding our numerical method for obtaining the same order of 
convergence for the partial derivatives of u as for u itself. In Section we present our 
scheme for the accurate computation of the Dirichlet to Neumann map. In Section 
the key elements of the algorithmic implementation of our method are discussed. We 
demonstrate the benehts of using our scheme in combination with an existing hnite 
element-based Grad-Shafranov solver in Section in which we focus on two magnetic 
conhnement fusion relevant examples. Finally, we summarize our work and suggest ideas 
for further development in Section 


2. Setup and General Idea 

The Grad-Shafranov equation is given by BE] 


d f 1 d'lp \ V’ 2^P 1 

dr \r dr) dz"^ 2 ’ 

with the components of the magnetic held in turn expressed as 

IdV’ H'd) IdV' 

r az r r or 


( 2 ) 

(3) 


In general, the right-hand side of Eq.([^ is a nonlinear function of "0, so that the equation 
has to be solved iteratively HD]. Without any penalty in terms of computational cost, 
one may simplify the left-hand side by the change of variables 




^/r' 


( 4 ) 
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This converts the Grad-Shafranov equation into the following semi-linear Poisson equa¬ 
tion: 


Am = 


dp{^/ru) 1 dP{^/ru) 3 u . 

- 3 - 7 ^ - 3 -+ 7 ^ •= ^ 


du 2r du 
where A denotes the cartesian-coordinate Laplacian - i.e. 


(5) 


dr'i Q^2 

and X := (r, z) G It is worth noting that the transformation Q is only advanta¬ 
geous in regions bounded away from the r = 0 axis. Fortunately, this requirement holds 
uniformly for domains corresponding to a large class of magnetic conhnement fusion 
applications, tokamaks in particular. 

We wish to solve equation ([^ in a smooth, bounded domain D, subject to homoge¬ 
neous Dirichlet boundary conditions - i.e. m = 0 on dD (which is equivalent to V’ = 0 on 
dD). To emphasize the generality of our approach, we de-emphasize the fusion specific 
variables and accordingly relabel the coordinates from (r, z) —>■ (x, y). A typical approach 
to solving ([^ is to use fixed-point iteration, computing by 


Am" = F(x,m^-i) , 


(7) 


with some specified initial guess vP. One stops the iteration when 


||Am"-F(x,m")|| <£, (8) 

for some reasonable norm and specified error tolerance £. When F(x, m) is proportional to 
M, a modihed version of the well-known inverse iteration method (TU] may be used to avoid 
the trivial solution m := 0 [iniini. Regardless of the iterative procedure, at each step, the 
Poisson equation ([^ has to be solved for m". While there are many techniques to solve ([^ 
the derivatives procedure introduced in this paper is independent of the specihc technique 
chosen. Stated differently our aim is to present a method for the accurate computation 
of derivatives that can be used with any existing hnite element-based Poisson or Grad- 
Shafranov code. Even so, in our numerical tests, we must choose a particular technique 
in order to obtain results. We use the isoparametric, bicubic Hermite finite element 
formulation used for the Grad-Shafranov equation in, for example, [SI El [TJ HH]. The 
method gives fourth order accuracy in u (and thus pj) and, as expected, third and second 
order accuracy in the first and second derivatives, respectively. We shall describe the 
chosen FEM scheme in sufficient detail in section 4 to understand the structure of the 
algorithm and the numerical results. We refer the reader interested in more details to 
the sources referenced above. 

The fundamental idea behind our technique is very simple. Say we wish to compute 
Ux- We differentiate (|^ to get 

Am^, - F„(x,M)Ma; = F„(x,m). (9) 


If the original numerical solution u is known, then this is a linear elliptic PDE for m^,, 
to which the same finite element method may be applied. In this way we in principle 
achieve the same order of accuracy as we did for m. This process may be repeated by 
further differentiation of the starting equation, at which point another linear elliptic 
PDE for the second derivative may be solved to the same order accuracy. Moreover, the 
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original equation (|^ was nonlinear, and thus required iteration over many FEM solves. In 
contrast, the equations for the derivatives are linear, requiring only one FEM solve. The 
computation of the derivatives is thus a small additional burden relative to the original 
computation. 

The one critical and computationally challenging point missing from the above dis¬ 
cussion is boundary data. In order to solve ([^, we need advance knowledge of on 
dD. Of course, knowing u to fourth order, we can easily use hnite differences to compute 
Ux on the boundary to third order, but this will lead to third order errors everywhere 
in the domain. We are back where we started - we know u to order k, but can only 
hnd Ux to order k — 1. However, there is a key difference now in that the function we 
need to differentiate is restricted to the boundary. Note that on the boundary, m, Mj., or 
any higher derivatives, are smooth, periodic functions of a single surface variable. It is 
thus an ideal target for spectral differentiation. If we can formulate the problem so that 
the only numerical derivatives required are spectral derivatives, we can hnd Ux on the 
boundary to the same order of accuracy as u. Admittedly, one may lose a constant scale 
factor in the accuracy of Ux relative to the accuracy of m, since spectral differentiation 
inherently increases the size of high frequency modes, but this is far preferable compared 
to a decrease in order of accuracy. 

In the following section, we demonstrate that through the use of integral equations 
and a clever change of variables, it is indeed possible to compute any derivatives of u 
on dD to the same order of accuracy as u itself while taking only spectral derivatives of 
periodic functions. 


3. Boundary Data for Derivatives 

To begin, it is convenient to decompose Vm on the boundary into components Un and 
Ut, which are normal and tangent to dD, respectively. Clearly, knowledge of the shape 
of dD allows Ux and Uy to be computed given Un and Ut- Throughout the article, we will 
use as our convention the inward normal direction and the counter-clockwise tangential 
direction unless otherwise noted. 

The tangential derivative is straightforward. Direct spectral differentiation of u re¬ 
stricted to the boundary yields Ut to the same order of accuracy as u. In our case, of 
course, = 0 is trivial to compute, but this will not be the case for second derivatives. 
The challenge lies in computing Intuitively, the difficulty comes from the fact that 
Un is sensitive to the behavior of u off the boundary, while Ut is not. This intuition is 
manifested mathematically by the fact that Green’s second identity gives 


/ G(x,x')m„(x') d/' = -m(x) 

'dD ^ 


(^(x, x')F(x, m(x')) dx' 


( 10 ) 


'D 


for any point x G dD, where G(x,x') = log(||x — x'||)/27r is the free space Green’s 
function of the Laplace operator. While this does give an integral equation for Un in 
terms of known quantities, as desired, it is an integral equation of the hrst kind - the 
unknown only appears inside an integral - which is ill-conditioned. Attempting to solve 
equation (10) for Un numerically is thus ill-advised at best, and in practice leads to very 


poor accuracy. 

What is needed is a transformation that converts normal derivatives to tangential 
derivatives, since the tangential derivatives are straightforward to compute spectrally. 
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For any harmonic function 0 (i.e. A0 = 0), there is just such a transformation: the 
harmonic conjugate of 0, which we denote by <F. It is dehned by 

= V0, (11) 

where V"*" = {—dy, dx)^■ Clearly, 0„ = <1)^, (pt = —and A0 = A<F = 0. The fact that 
0 is harmonic is essential to this transformation; if A0 ^ 0, then the mixed partials of <I> 
cannot be equal, indicating the nonexistence of such a 

Unfortunately, u is not harmonic. However, it is easily decomposed into harmonic 
and anharmonic parts. We write u = where is given by 


M^(x) = / G(x, x')F(x', m(x')) dx', (12) 

J D 

Here again G(x,x') = log(||x — x'||)/27r is the free-space Green’s function of the Laplace 
operator. We often abbreviate this G when the arguments can be understood. The 
function is the anharmonic contribution to the total u. Meanwhile, the harmonic 
contribution corresponds to the homogeneous solution and satishes 


Au^ = 0 , 


u 


IdD 


= 


IdD ■ 


(13) 


We can now separate the problem into the computation and Importantly, is 
harmonic, so it has a harmonic conjugate satisfying = u^. The strategy is thus to 
compute on dD and to take its spectral tangential derivative to hnd u^. To hnd wjj, 
we differentiate (12) analytically and evaluate the resulting integral numerically. Note 
hnally that Eq. (13) is relatively simple because u has homogeneous boundary data. In 
general, the boundary condition on would be u^\qj^ = '^\dD ~ '^^\dD- general 

form will be necessary when computing higher order derivatives. 

The following two subsections lay out the details of computing wjj and Since 
evaluating wf* - which can be done concurrently with evaluation of - turns out to be 
a key ingredient in Ending «(), we first describe the procedure for wjj first. In all that 
follows, the goal is to find values of wjj and at points Xj G dD that are evenly spaced 
in the arc-length variable s. The details of computing the arc-length grid are discussed 
in appendix A. 


3.1. Computing 

As just mentioned, we will eventually need as well as wjj. We evaluate both from 
the exact formula 

[ VGF{x',u{x'))dx', (14) 

J D 


where 


VG(x,x') 


1 x-x' 
27r ||x — x'l 


(15) 


A major advantage of the Green’s function formulation is that we have an exact inte¬ 
gral representation for the partial derivatives of given by Eq. This formulation. 


however, comes with two well known computational challenges: first, a robust high order 
quadrature rule must be used to evaluate the singular integrals; second, the quadrature 
can be potentially costly in terms of computational time since the integrals are over a 
two-dimensional domain. Our numerical method addresses both challenges by relying on 
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D 


Figure 1: Domain D and its exterior E. c is an expansion center corresponding to xq in the 
exterior. 


a variant of the Quadrature By Expansion (QBX) quadrature scheme for the singular 
integrals [12] accelerated by the Fast Multipole Method (EMM) [2U] . 

Specihcally, let us £x our attention to evaluating Vup at xq G dD. Let E = M.‘^\D he 
the exterior of D. We observe that both the components of Vup are harmonic in E and 
hence dehne smooth functions there. Using a smooth quadrature rule, we can evaluate 
Vwp and its derivatives accurately at c = xq + rno, with r a constant such that r = 0{h) 
and h is the local spacing of discretization points on the boundary (see Figure [^. To 
make this more precise, let x = (x, y), x.' = (x', y'), c = (c^,, Cy), z = x + iy, z' = x' + iy' 
and c = Cx -\- iCy. We observe that w = dxU^ — idyU^ is complex analytic in E, since 
the real and imaginary parts are harmonic and complex conjugates of each other. Using 


equation (14), we get. 


w = 


E{z', u{z')) 


dx . 


(16) 


2n Jz — z' 

Since w is complex analytic in E, we can form a Taylor expansion about c, to obtain the 
following representation for to. 


w = ^{z 

j=0 


E{z', u{z')) 


dx . 


Id (U — cY 

It can be shown that a pth order truncated expansion of the above equation, given by 

F(z'Mz')) 


( 17 ) 


U -Id (z’-cy 


is a high order approximation to w, even at xq. More precisely, 

\w (xo) - w (X(i)l < C (p, D) r”. 


(18) 


(19) 


We refer the reader to nain] for a detailed discussion. To evaluate w, we still need to 
compute the integrals, 

( F(z'Mz')) ^ ( 20 ) 

Id 


(z' — cY 


which can be done to high order using a smooth quadrature rule, since the integrand is 
now smooth. We choose 4th order tensor-product Gauss-Legendre quadrature on each 
interior grid cell - and 8th order on each boundary cell - in mapped polar-like coordinates 


(see (30)). 

The computational cost of a naive implementation of the above scheme would be 
O {Nh ■ Nx), where W is the number of points on the boundary where we wish to evaluate 
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and is the number of volume points used to discretize the integrands in equation 
(20). However, the above computation can be accelerated by the FMM as discussed 
in [221 E3| to reduce the computational cost to O {Nb + N.^), which is of comparable 
complexity to the other steps in our method. 

3.2. Computing 

Recall that to evaluate our approach is to compute U^, the harmonic conjugate 
of u^. Since is harmonic by construction, applying Green’s second identity to it gives 


tf/^x) = [ (G!/*(x') - G„(/'‘(x')) dl'. 

^ JdD 


IdD 

Th _ ^,h _ r,,P 


By rearranging and noting that = —uf on the boundary, we have 

1 , 


-HRx)+ / GnU^U)dl' = 


'an 


I an 


Gu^ (x') dl' 


( 21 ) 


( 22 ) 


for any x G dD. 

Given the computation of according to (14), it is easy to evaluate uf, so that 
we can regard (22) as a second-kind integral equation for the unknown We reiterate 
that it is crucial that this integral equation is of the second kind, in contrast to (10). 
This integral equation can thus be accurately solved by discretizing each of the integrals 
on our arc-length grid. It is worth mentioning that the integrand on the left-hand side of 
(22) is not singular in spite of its appearance. This is because 


lim G„(x,x') = —k(x), 
47r 


(23) 


where k is the signed curvature of dD. Given a parameterization of the curve {x{t), y{t)), 
it is dehned by 

P4) 

Thus, Gn restricted to the boundary is a smooth function once we dehne Gn(x,x) := 
K(x)/47r. As such, we are free to dehne 

Pij . As Gjj(xj, Xj), (25) 

where As is the arc-length distance between grid points. Thus, the approximation of the 
integral on the left-hand side of (22) by the trapezoidal rule is written as 


[ GnU\x')di' 
JaD 


(26) 


The trapezoidal rule suffices since it is well known to converge spectrally for smooth, 
periodic functions m- 

The integrand on the right-hand side of (22), on the other hand, is logarithmically 
singular. We denote by the value of the integral at x = Xj, and again compute it using 
the standard version of QBX described above. 

In the end, the discretized version of (22) is 

E ( +«<i) 


( 27 ) 











This linear system is ill-conditioned because is defined in terms of its derivatives, and 


thus only well defined up to a constant. We thus expect infinitely many solutions to (27). 
The problem is solved by imposing the additional constraint that average to zero. 
That is, 




= 0 . 


(28) 


We may impose this constraint without increasing the size of the system by solving 


E 


1 

2 Op- + Pij 
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(29) 


As shown in [25], this new linear system has a unique solutions which also solves Eq. (22) 
with probability 1. 


The linear system (29) may now be solved using any standard linear algebra package. 


It bears mentioning that this system is dense, while the stiffness matrix in the FEM 
formulation is sparse. It may thus seem that solving (29) becomes the rate limiting 
step in the computation. However, this is not the case. The dimension of the dense 
system scales with the number of grid points on the boundary, while the stiffness matrix 
size scales with the number of points in the entire two dimensional domain - roughly the 
square of the number of boundary points. Thus, the complexity of solving each of the two 


linear systems is comparable. Moreover, the system (29) need only be solved once, while 


the finite element system must in general be solved several times since the semi-linear 
Poisson equation ([^ has to be solved iteratively. 


4. Implementation Details and Algorithm Summary 

As mentioned in Section we choose an isoparametric, bicubic Hermite finite element 
formulation, which we briefly outline here. The computational domain is represented in 
(p, 6) coordinates, which are defined in terms of the Cartesian coordinates x and y by the 
following transformation 


X = Xc + pf{0) cos 9 

y = Vc + pfW sin 9 


(30) 


for some pre-specified central axis {xc,yc)- Accordingly, 9 is the polar angle, and f{9) 
is represented by a Fourier series, which is computed from given boundary information. 
When necessary, tangent and normal directions to the boundary may be found by differ¬ 
entiating this Fourier series. 

In {p,9) coordinates, the computational domain is [0,1] x [0,27r). This is subdivided 
into an N X N grid, to which the FEM is applied using the bicubic Hermite basis in the 
(p, 9) variables. The integrals of the basis functions that must be evaluated are computed 
using tensor-product Gauss-Legendre quadrature with four nodes in each direction. The 
solution M of (|^ is approximated by iteratively applying this FEM scheme to ([^, stopping 
according to the criterion ([^, using infinity norm over the unknowns and e = 1.7 x 10“^^. 

For the purposes of spectral differentiation, it is desirable to have a grid along the 
boundary that is evenly spaced in arc length. The angles 9j corresponding to such a 
grid may be computed using the method described in appendix A. We construct this 
grid using 8N points, and denote the arc-length grid-points on the boundary by x^. We 
proceed through the following steps to find and Uy on the boundary. Throughout, 
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rij = denotes the unit inward normal to the boundary at xj, and tj the unit 

counterclockwise tangent vector. We present the method in the more complex case in 
which it is not assumed that m = 0 on the boundary, for this general case is required 
when computing second derivatives. 

1. Use FMM accelerated QBX to evaluate 


Vmj = ^ VG(xj,x')F(x',m(x')) dx', (31) 

Jd 

as described in 3.1. 

2. Set = iij ■ and = tj ■ Vm^(xj). 

3. Compute 

Jj = [ G(xj,x') [Mi(x') - M?(x')] d/' (32) 

JdD 

by again using FMM accelerated QBX applied to the computed values of u^j. Here, 
Ut is computed by direct spectral differentiation of the FEM solution. 

4. Solve the linear system 

Uj = 'ji, (33) 

where As is the spacing of the arc-length grid, as described in 3.2. In our code, we 
simply use MATLAB’s “backlash” operator, as we hnd this is not the rate limiting 
step in the procedure. More generally, an iterative solver based on the generalized 
minimal residual method (GMRES) is satisfying since the linear system only needs 
to be inverted once. The system is well conditioned so will converge quickly, i.e in 
0(1) iterations. 

5. Compute by spectral differentiation of U^. 

6. Set Unj = U^j + u^j, from which it follows that = n^jUn,] + UyjUtj, and 

'^yd ~ ^yd'^ri,j ~ 'nxjUtj. 

With this boundary data in hand, we use the same hnite element formulation to solve 


Aux - E«(x, u)ux = Fxi^, u) 
AUy - F„(x,M)My = iS/(x,M) 


(34) 


for Ux and Uy, respectively. 

Computing second derivatives is directly analogous. The procedure above is re-used. 


but with u ^ Ux everywhere and F(x, u) —)■ Fx{'x.,u) -|- Fu{'k,u)ux in (31). To be more 
specihc, we evaluate 


V< = / VG(x,x') [Fa.(x',M(x')-h F„(x',m(x ))m3.(x')] dx 


(35) 


’ D 


on the arc-length grid, and solve 


[ GnU^{x)dl'= [ G(x,x') [Mxt(x')-<t(x')] d/' (36) 

^ JdD JdD 

for U!^. Then, Uxn = + Uxn is found by direct spectral differentiation of Ux on 

the same grid. 


10 



Having u^n and Uxt-, we may find u^x and u^y on the boundary, which are used as 
boundary data to solve 


^)^XX Fxx H“ ‘^Fxu^X F Fxu^x 
^'^xy Fxi{^^ U^Uxy Fxy H“ Fxu^y H“ FyyVjx FxyUxUy 


(37) 


for Uxx and Uxy respectively. It is then simple to find Uyy, since 


'^yy Uxx 

as a consequence of the original PDE (|^. Once u and its derivatives are known, it is 
straightforward to compute V’ and its derivatives to the same accuracy using the formula 

0. 


5. Numerical Results 


For numerical tests, we consider two situations in which exact solutions are known. 
The first situation corresponds to toroidally axisymmetric plasma equilibria with pressure 
and current profiles chosen so that the right-hand side of the Grad-Shafranov equation 
does not depend on ^fJ. Exact solutions are compared with numerical solutions for two 
magnetic conhnement devices, as discussed in detail in section 5.1 In this hrst situation. 


both equations ([^ and ([^ are linear partial differential equations. In order to show 
that our procedure applies just as well to nonlinear partial differential equations, we 


consider in section 5.2 an exact solution to the nonlinear Poisson-Boltzmann equation on 


an ellipse-like domain. 


5.1. Grad-Shafranov equation with Solov’ev profiles 

In this hrst example, we solve the Grad-Shafranov equation with a Solov’ev pressure 
prohle [2ni EZI and no diamagnetic or paramagnetic contribution to the toroidal magnetic 
held {dP/dip = 0). The simplihed Grad-Shafranov equation is given by 


d f 1 dfi\ d'^ip 
dx \x dx J dy'^ 


Cx^ 


(39) 


A simple and exact solution to this equation is 




C 

J 


x^ -\- di-\- d2x‘^ -t- dfix"^ 


Ax^y"^) 


(40) 


for any constants di, d 2 , and ds, where the boundary dD corresponding to "0 = 0 is 
simply chosen to be the curve along which the above polynomial vanishes. Of course, not 
every choice of G, di, d 2 , and da results in a reasonable plasma boundary. We choose 
these constants following the approach in uni [23. Namely, we rewrite them in terms of 
plasma relevant quantities. These are e, d, and k, which are called the inverse aspect 
ratio, triangulation, and elongation, respectively. By imposing 


+ e, 0) = — e) = fiil — de, ne) = 0, (41) 

we may write an invertible linear system that relates (e, d, k) to {di,d 2 ,d^) [10]. In all 
our examples, we will use C = 10. Its value is not important since it can be scaled out 
of the problem. 
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We compute the errors in u and its derivatives using the norm. That is, the error 
in any quantity Q(x) is given by 


Error 



(42) 


The integral is evaluated using Gauss-Legendre quadrature over the FEM grid in (s, 0) 
coordinates. In addition, we measure error in the location of the magnetic axis, which is 
the value of x at which Vi/^ = 0. This is computed via gradient descent on the y = 0 axis, 
since symmetry tells us the magnetic axis must lie on this line. The numerical result is 
compared to the exact axis location, given by 


X 


* 


-do 


C + 8^3 


(43) 


We consider two cases: 1) (e, 5, k) = (0.32,0.33,1.7), which corresponds to the ge¬ 
ometry of the ITER tokamak [28]; 2) (e, 5, k) = (0.78,0.35,2), which corresponds to the 
geometry of the spherical tokamak NSTX [22] • For reference, the equilibria are shown in 
hgure Since the equations for the Solov’ev equilibria are in fact linear, we perform an 
additional test on a Poisson-Boltzmann equation to conhrm the method’s performance 
on nonlinear equations. 



Figure 2: Left: Contours of if for the ITER-like Solov’ev equilibrium given by Eq.(39) with 
(e,(5, k) = (0.32,0.33,1.7). Right: Contours of if for the NSTX-like Solov’ev equilibrium given 
by Eq.(39) with {e,S,K) = (0.78,0.35,2). 


5.1.1. ITER example 

Let us hrst focus on the ITER-like equilibrium. In hgure we plot the error in the 
solution u. As expected, we see that this error is 0{N~^) as a result of the FEM scheme 
used. More interestingly, in hgure we plot the error in and Uy, as well as the location 
of the magnetic axis. The dashed lines are computed by reading oh uq and Ug from the 
hnite element solution, and then converting to rectangular coordinates. As expected, one 
order of accuracy is lost relative to the solution. 

In contrast, the solid lines are the errors using the new method presented in this 
article. One can observe the improved convergence rate and improved absolute error for 
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Solution Error 



Figure 3: 
(e, 5, k) = 


Error in the solution u for the ITER-like Solov’ev equilibrium given by Eq.(39) with 
(0.32,0.33,1.7), displaying the expected N~^ convergence rate. 


1 St Derivative Error 




Figure 4- Left: Error in the first partial derivatives for the ITER-like Solov’ev equilibrium given 
by Eq.(39) with (e, <5, k) = (0.32, 0.33,1.7). Dashed lines represent the naive computation, while 
results obtained with the new method presented in this article are shown in solid lines. Right: 
Error in the location of the magnetic axis, displaying improved convergence relative to results 
directly obtained with CREASE and reported in mn\. 


any N > 32. At iV = 256, we see an improvement of more than an order of magnitude in 
each derivative. The error in the magnetic axis is also 0{N~'^). This is to be compared 
with the results presented in [6] , which did not take advantage of our method, and where 
0{N~^) accuracy was observed using the same finite element basis. 

Even more signihcant are the improvements in accuracy of the second derivatives, 
shown in figure For the clarity of the figure, we did not plot the error in Uyy obtained 
with the new method because it is similar to the error in Uxx according to equation 
(38). We did not plot the error in Uxy that one obtains with the standard finite element 
approximation either because we could not get it to converge properly. We note by 
looking at figure that the improved convergence rate with the new method is observed 
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2 nd Derivative Error 



Figure 5: Error in the second derivatives for the ITER-like Solov’ev equilibrium given by Eq.(39) 
with (e,5, k) = (0.32,0.33,1.7). Dashed lines represent the naive computation, while results 
obtained with the new method presented in this article are shown in solid lines. 


for moderate iV, but the curve flattens as N grows beyond ~ 128. This is due to round-off 
error that accumulates from taking numerous spectral derivatives. Even so, at the present 
convergence rate, the standard hnite element approximation requires N 5000 to reach 
the accuracy the new method achieves at iV = 256. That is a reduction in complexity 
of roughly (5000/256)^ ~ 381, which easily compensates for the additional computation 
resulting from execution of the methods presented here. An alternate comparison is that 
for N = 128 the new method is more than 1000 times more accurate for u^x than the 
direct FEM result. 

Even in the event that only 4 digits of accuracy are required in the second derivatives, 
the new scheme accomplishes this with 50, while the standard method requires 

N K. 315. The speed improvement in this case is roughly (315/50)^ ~ 40. This again 
easily outstrips the additional cost devoted to the boundary and the extra hnite element 
solves. 


5.1.2. NSTX example 

As a second numerical example, we consider the more challenging case of a Solov’ev 
equilibrium as given by Eq. ( [M| ) for the high inverse aspect ratio, highly elongated spher¬ 
ical tokamak NSTX [22], with typical parameters (e, 5, k) = (0.78,0.35,2). The results 


are shown in hgures and We did not plot the curves corresponding to Uyy with the 
new method and Uxy with the standard method for the same reasons as the ones we gave 


in section 5.1.1 We observe the same improvement in terms of the order of convergence 


of the derivatives of the solution to the Grad-Shafranov equation. One does note that for 
the magnetic axis and second derivatives, the convergence is not as smooth at small N as 
in the ITER case. We hypothesize that this is due to the fact that the NSTX boundary is 
more difficult to resolve so that larger N is required to observe the asymptotic behavior. 
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Figure 6: Left: Error in the first partial derivatives for the NSTX-like Solov’ev equilibrium 
given by Eq.(39) with (e, <5, k) = (0.78,0.35,2). Dashed lines represent the naive computation, 
while results obtained with the new method presented in this article are shown in solid lines. 
Right: Error in the location of the magnetic axis. 


2nd Derivative Error 



Figure 1: Error in the second derivatives for the NSTX-like Solov’ev equilibrium given by 
Eq.(39) with (e, <5, k) = (0.78,0.35,2). Dashed lines represent the naive computation, while 
results obtained with the new method presented in this article are shown in solid lines. 


5.2. Poisson-Boltzmann example 

We now verify that the method we describe in this article also performs well for 
nonlinear cases by considering the nonlinear Poisson problem given by 



Au = ae '^IdD — 0) 

(44) 

where dD is dehned by 

Cl cosh ky — C 2 cos kx = 1. 

(45) 

We pick a = 2kf{c\ — c\), 

giving an exact solution, 



u = 2 log(ci cosh ky — C 2 cos kx), 

(46) 
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against which our numerical solutions may be compared. We evaluate the numerical error 
in the same way as explained in section 5.1 For the free parameters, we choose fc = 7r/5, 
Cl = .0287 and C 2 = .3301. This gives an ellipse-like boundary with a 2-to-l aspect ratio. 

Results for the hrst and second derivatives are plotted in hgure The method 
presented here performs as expected. We again observe the flattening of the second 
derivative curve due to accumulated round-off error from numerous spectral derivatives. 
Curiously, the hnite difference approximation of second derivatives behaves even worse 
than expected. We’ve plotted Uxx, which initially converges at the expected rate before 
diverging at large N. However, Uxy and Uyy fail to converge at all when hnite differences 
are used. This hints that in addition to being more accurate, the new method may prove 
more robust than the standard approach. 


1 St Derivative Error 2nd Derivative Error 




Figure 8: Left: Error in the first partial derivatives for the Poisson-Boltzmann example. Dashed 
lines represent the naive computation, while results obtained with the new method presented 
in this article are shown in solid lines. Right: Error in the second derivatives, with dashed lines 
again corresponding to the naive computation and solid lines to the new method. 


6. Discussion and Conclnsions 

We have presented a method based on an integral equation formulation for the accu¬ 
rate numerical evaluation of the Dirichlet to Neumann map for Poisson’s equation and 
the Grad-Shafranov equation. The method takes as an input the solution to either of 
these partial differential equations obtained with any hnite element-based solver. By dif¬ 
ferentiating the partial differential equation analytically, we are then able to use the same 
hnite element solver to compute partial derivatives of the solution with the same order of 
convergence as the solution itself. The computational cost of implementing the numer¬ 
ical method we describe in this article is comparable to two Poisson or Grad-Shafranov 
solves. This additional computational cost is ohset by the fact that on small grids, we 
obtain accuracies for the hrst and second derivatives that could only be achieved on a 
much larger grid with the standard method. In the context of plasma physics simulations, 
a code based on our numerical method can be advantageously combined with existing 
and commonly used Grad-Shafranov solvers to calculate quantities such as the magnetic 
held, the parallel current density and the magnetic curvature with much higher accuracy 
than is possible with the standard methods on the fairly small grids that can usually 
be ahorded. At present, our method is only applicable to cases for which the boundary 
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of the domain is smooth. The quadrature schemes we employ here can handle bound¬ 
aries with a sharp corner, but the Fourier based method we use to compute spectrally 
accurate derivatives on the boundary would not yield accurate results for curves that are 
not smooth. Given the importance of boundaries with corners in magnetic conhnement 
fusion, an improved scheme with this capability, most likely based on local interpolation 
instead of global interpolation, would be desirable. This is the subject of ongoing work, 
with progress to be reported at a later date. 
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Appendix A. Arc-length grid computation 

We desire a sequence of angles 9j that correspond to equally spaced points in arc- 
length, with spacing given by As = 8N/L, where L is the total length of the curve. We 
first compute L by evaluating 


L = 




(A.l) 


where dx/dO and dy/dO are computed by differentiating (30), using the Fourier series 
representation of / to compute its derivative. This integral can be evaluated accurately 
using any high order quadrature desired. We use Gauss-Legendre with 1000 nodes. 

One then sets = 0, and solves 


8N 

~T 




(A.2) 


for 9j given 9j-i using Newton’s method. The integral is again evaluated using any 
desired quadrature. We use Gauss-Legendre with 16 nodes. 
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